#include "SymmetryBoundary.h"
#include "FieldCircuitCoupling.h"
#include <math.h>
#include <stdexcept>
#include <algorithm>
#include <map>

namespace SymmetryBoundary {

SymmetryBoundaryManager::SymmetryBoundaryManager() {
    // Initialize symmetry planes with default values (disabled)
    symmetryPlanes.resize(3);
    
    // XY plane (z=0)
    symmetryPlanes[0].plane = PlaneOrientation::XY_PLANE;
    symmetryPlanes[0].type = SymmetryType::NONE;
    symmetryPlanes[0].planeOffset = 0.0;
    symmetryPlanes[0].enabled = false;
    
    // YZ plane (x=0)
    symmetryPlanes[1].plane = PlaneOrientation::YZ_PLANE;
    symmetryPlanes[1].type = SymmetryType::NONE;
    symmetryPlanes[1].planeOffset = 0.0;
    symmetryPlanes[1].enabled = false;
    
    // XZ plane (y=0)
    symmetryPlanes[2].plane = PlaneOrientation::XZ_PLANE;
    symmetryPlanes[2].type = SymmetryType::NONE;
    symmetryPlanes[2].planeOffset = 0.0;
    symmetryPlanes[2].enabled = false;
}

SymmetryBoundaryManager::~SymmetryBoundaryManager() {
}

void SymmetryBoundaryManager::setSymmetryConfig(const SymmetryConfig& config) {
    int index = -1;
    switch(config.plane) {
        case PlaneOrientation::XY_PLANE:
            index = 0;
            break;
        case PlaneOrientation::YZ_PLANE:
            index = 1;
            break;
        case PlaneOrientation::XZ_PLANE:
            index = 2;
            break;
        default:
            throw std::invalid_argument("Invalid plane orientation");
    }
    
    symmetryPlanes[index] = config;
}

bool SymmetryBoundaryManager::getSymmetryConfig(PlaneOrientation plane, SymmetryConfig& config) const {
    for (const auto& symmetryPlane : symmetryPlanes) {
        if (symmetryPlane.plane == plane) {
            config = symmetryPlane;
            return true;
        }
    }
    return false;
}

bool SymmetryBoundaryManager::applySymmetryToImpedanceMatrix(
    std::vector<std::vector<FieldCircuitCoupling::Complex>>& impedanceMatrix,
    const std::vector<Vector_3D>& basisFunctionPositions) const {
    try {
        // Apply symmetry for each enabled plane
        for (const auto& config : symmetryPlanes) {
            if (config.enabled && config.type != SymmetryType::NONE) {
                if (!applySinglePlaneSymmetry(impedanceMatrix, basisFunctionPositions, config)) {
                    return false;
                }
            }
        }
        
        return true;
    } catch (...) {
        return false;
    }
}

bool SymmetryBoundaryManager::applySymmetryToExcitation(
    std::vector<FieldCircuitCoupling::Complex>& excitation,
    const std::vector<Vector_3D>& basisFunctionPositions) const {
    try {
        // Apply symmetry transformations to the excitation vector
        std::vector<FieldCircuitCoupling::Complex> tempExcitation = excitation;
        
        for (const auto& config : symmetryPlanes) {
            if (config.enabled && config.type != SymmetryType::NONE) {
                for (size_t i = 0; i < basisFunctionPositions.size(); ++i) {
                    // Check if the basis function is on the positive side of the symmetry plane
                    double coord = getPlaneCoordinate(basisFunctionPositions[i], config.plane);
                    
                    if (coord > config.planeOffset + 1e-9) {  // Avoid floating point errors
                        // Find the image basis function
                        Vector_3D imagePos = computeImagePosition(basisFunctionPositions[i], config);
                        
                        // Find the index of the image basis function
                        for (size_t j = 0; j < basisFunctionPositions.size(); ++j) {
                            // Check if j is the image of i
                            if (fabs(basisFunctionPositions[j].x - imagePos.x) < 1e-9 &&
                                fabs(basisFunctionPositions[j].y - imagePos.y) < 1e-9 &&
                                fabs(basisFunctionPositions[j].z - imagePos.z) < 1e-9) {
                                
                                // Apply symmetry factor to the image excitation
                                Vector_3D dirVector = {1.0, 1.0, 1.0};  // Default direction
                                FieldCircuitCoupling::Complex factor = getSymmetryFactor(config.type, dirVector);
                                
                                // For certain symmetry types, we need to combine excitations
                                if (config.type == SymmetryType::SYMMETRIC) {
                                    tempExcitation[i] = (excitation[i] + factor * excitation[j]) * 0.5;
                                } else if (config.type == SymmetryType::ANTISYMMETRIC) {
                                    tempExcitation[i] = (excitation[i] - factor * excitation[j]) * 0.5;
                                }
                                break;
                            }
                        }
                    }
                }
            }
        }
        
        excitation = tempExcitation;
        return true;
    } catch (...) {
        return false;
    }
}

bool SymmetryBoundaryManager::expandSolution(
    std::vector<FieldCircuitCoupling::Complex>& solution,
    const std::vector<Vector_3D>& basisFunctionPositions) const {
    try {
        // Create a copy of the solution to work with
        std::vector<FieldCircuitCoupling::Complex> expandedSolution = solution;
        
        // For each symmetry plane, expand the solution by considering image currents
        for (const auto& config : symmetryPlanes) {
            if (config.enabled && config.type != SymmetryType::NONE) {
                for (size_t i = 0; i < basisFunctionPositions.size(); ++i) {
                    double coord = getPlaneCoordinate(basisFunctionPositions[i], config.plane);
                    
                    if (coord > config.planeOffset + 1e-9) {  // Avoid floating point errors
                        // Find the image basis function
                        Vector_3D imagePos = computeImagePosition(basisFunctionPositions[i], config);
                        
                        for (size_t j = 0; j < basisFunctionPositions.size(); ++j) {
                            if (fabs(basisFunctionPositions[j].x - imagePos.x) < 1e-9 &&
                                fabs(basisFunctionPositions[j].y - imagePos.y) < 1e-9 &&
                                fabs(basisFunctionPositions[j].z - imagePos.z) < 1e-9) {
                                
                                // Determine the direction vector for the basis function
                                Vector_3D dirVector = {1.0, 1.0, 1.0};  // Default
                                
                                // Apply symmetry factor to determine the image current
                                FieldCircuitCoupling::Complex factor = getSymmetryFactor(config.type, dirVector);
                                
                                // Set the image current based on symmetry type
                                if (config.type == SymmetryType::SYMMETRIC) {
                                    expandedSolution[j] = factor * solution[i];
                                } else if (config.type == SymmetryType::ANTISYMMETRIC) {
                                    expandedSolution[j] = -factor * solution[i];
                                } else if (config.type == SymmetryType::PEC) {
                                    // For PEC, tangential electric field is zero
                                    expandedSolution[j] = solution[i];
                                } else if (config.type == SymmetryType::PMC) {
                                    // For PMC, tangential magnetic field is zero
                                    expandedSolution[j] = -solution[i];
                                }
                                break;
                            }
                        }
                    }
                }
            }
        }
        
        solution = expandedSolution;
        return true;
    } catch (...) {
        return false;
    }
}

Vector_3D SymmetryBoundaryManager::computeImagePosition(
    const Vector_3D& point, 
    const SymmetryConfig& config) const {
    Vector_3D imagePoint = point;
    double coord = getPlaneCoordinate(point, config.plane);
    double mirroredCoord = 2 * config.planeOffset - coord;
    
    return setPlaneCoordinate(point, config.plane, mirroredCoord);
}

bool SymmetryBoundaryManager::isAnySymmetryEnabled() const {
    for (const auto& config : symmetryPlanes) {
        if (config.enabled && config.type != SymmetryType::NONE) {
            return true;
        }
    }
    return false;
}

FieldCircuitCoupling::Complex SymmetryBoundaryManager::getSymmetryFactor(
    SymmetryType type,
    const Vector_3D& vector) const {
    switch (type) {
        case SymmetryType::SYMMETRIC:
        case SymmetryType::PEC:
            return FieldCircuitCoupling::Complex(1.0, 0.0);
        case SymmetryType::ANTISYMMETRIC:
        case SymmetryType::PMC:
            return FieldCircuitCoupling::Complex(-1.0, 0.0);
        default:
            return FieldCircuitCoupling::Complex(1.0, 0.0);
    }
}

bool SymmetryBoundaryManager::initialize(std::shared_ptr<FieldCircuitCoupling::FieldSolverInterface> solver) {
    try {
        fieldSolver = solver;
        return true;
    } catch (...) {
        return false;
    }
}

int SymmetryBoundaryManager::getReductionFactor() const {
    int factor = 1;
    
    // For each enabled symmetry plane, multiply the reduction factor
    for (const auto& config : symmetryPlanes) {
        if (config.enabled && config.type != SymmetryType::NONE) {
            factor *= 2;  // Each symmetry plane reduces problem size by half
        }
    }
    
    return factor;
}

bool SymmetryBoundaryManager::generatePositionMap(
    const std::vector<Vector_3D>& originalPositions,
    std::vector<std::pair<size_t, size_t>>& positionMap) const {
    try {
        positionMap.clear();
        
        // For each enabled symmetry plane, create a mapping between original and image positions
        for (const auto& config : symmetryPlanes) {
            if (config.enabled && config.type != SymmetryType::NONE) {
                for (size_t i = 0; i < originalPositions.size(); ++i) {
                    double coord = getPlaneCoordinate(originalPositions[i], config.plane);
                    
                    if (coord > config.planeOffset + 1e-9) {  // Positive side of the plane
                        Vector_3D imagePos = computeImagePosition(originalPositions[i], config);
                        
                        // Find the corresponding image position index
                        for (size_t j = 0; j < originalPositions.size(); ++j) {
                            if (i != j &&
                                fabs(originalPositions[j].x - imagePos.x) < 1e-9 &&
                                fabs(originalPositions[j].y - imagePos.y) < 1e-9 &&
                                fabs(originalPositions[j].z - imagePos.z) < 1e-9) {
                                
                                positionMap.push_back(std::make_pair(i, j));
                                break;
                            }
                        }
                    }
                }
            }
        }
        
        return true;
    } catch (...) {
        return false;
    }
}

// Private methods
bool SymmetryBoundaryManager::isPointOnSymmetryPlane(
    const Vector_3D& point, 
    const SymmetryConfig& config) const {
    double coord = getPlaneCoordinate(point, config.plane);
    return fabs(coord - config.planeOffset) < 1e-9;
}

double SymmetryBoundaryManager::getPlaneCoordinate(
    const Vector_3D& point, 
    PlaneOrientation plane) const {
    switch (plane) {
        case PlaneOrientation::XY_PLANE:
            return point.z;
        case PlaneOrientation::YZ_PLANE:
            return point.x;
        case PlaneOrientation::XZ_PLANE:
            return point.y;
        default:
            throw std::invalid_argument("Invalid plane orientation");
    }
}

Vector_3D SymmetryBoundaryManager::setPlaneCoordinate(
    const Vector_3D& point, 
    PlaneOrientation plane, 
    double value) const {
    Vector_3D result = point;
    
    switch (plane) {
        case PlaneOrientation::XY_PLANE:
            result.z = value;
            break;
        case PlaneOrientation::YZ_PLANE:
            result.x = value;
            break;
        case PlaneOrientation::XZ_PLANE:
            result.y = value;
            break;
        default:
            throw std::invalid_argument("Invalid plane orientation");
    }
    
    return result;
}

bool SymmetryBoundaryManager::applySinglePlaneSymmetry(
    std::vector<std::vector<FieldCircuitCoupling::Complex>>& impedanceMatrix,
    const std::vector<Vector_3D>& basisFunctionPositions,
    const SymmetryConfig& config) const {
    try {
        size_t n = impedanceMatrix.size();
        
        // Create a temporary matrix to store the modified impedance
        std::vector<std::vector<FieldCircuitCoupling::Complex>> tempMatrix = impedanceMatrix;
        
        for (size_t i = 0; i < n; ++i) {
            for (size_t j = 0; j < n; ++j) {
                // Check if both basis functions are on the positive side of the plane
                double coordI = getPlaneCoordinate(basisFunctionPositions[i], config.plane);
                double coordJ = getPlaneCoordinate(basisFunctionPositions[j], config.plane);
                
                // Skip if either is on the plane or in the negative side
                if (coordI <= config.planeOffset + 1e-9 || coordJ <= config.planeOffset + 1e-9) {
                    continue;
                }
                
                // Compute image positions
                Vector_3D imagePosI = computeImagePosition(basisFunctionPositions[i], config);
                Vector_3D imagePosJ = computeImagePosition(basisFunctionPositions[j], config);
                
                // Find indices of image basis functions
                size_t imageI = n, imageJ = n;
                for (size_t k = 0; k < n; ++k) {
                    if (fabs(basisFunctionPositions[k].x - imagePosI.x) < 1e-9 &&
                        fabs(basisFunctionPositions[k].y - imagePosI.y) < 1e-9 &&
                        fabs(basisFunctionPositions[k].z - imagePosI.z) < 1e-9) {
                        imageI = k;
                    }
                    if (fabs(basisFunctionPositions[k].x - imagePosJ.x) < 1e-9 &&
                        fabs(basisFunctionPositions[k].y - imagePosJ.y) < 1e-9 &&
                        fabs(basisFunctionPositions[k].z - imagePosJ.z) < 1e-9) {
                        imageJ = k;
                    }
                }
                
                // If both image indices are valid
                if (imageI != n && imageJ != n) {
                    // Determine symmetry factors based on boundary type
                    FieldCircuitCoupling::Complex factorII, factorJJ, factorIJ, factorJI;
                    
                    switch (config.type) {
                        case SymmetryType::SYMMETRIC:
                            factorII = factorJJ = factorIJ = factorJI = FieldCircuitCoupling::Complex(1.0, 0.0);
                            break;
                        case SymmetryType::ANTISYMMETRIC:
                            factorII = factorJJ = FieldCircuitCoupling::Complex(1.0, 0.0);
                            factorIJ = factorJI = FieldCircuitCoupling::Complex(-1.0, 0.0);
                            break;
                        case SymmetryType::PEC:
                            // For PEC, tangential E field is zero
                            factorII = factorJJ = FieldCircuitCoupling::Complex(1.0, 0.0);
                            factorIJ = factorJI = FieldCircuitCoupling::Complex(-1.0, 0.0);
                            break;
                        case SymmetryType::PMC:
                            // For PMC, tangential H field is zero
                            factorII = factorJJ = FieldCircuitCoupling::Complex(1.0, 0.0);
                            factorIJ = factorJI = FieldCircuitCoupling::Complex(1.0, 0.0);
                            break;
                        default:
                            continue;
                    }
                    
                    // Apply symmetry to the impedance matrix
                    tempMatrix[i][j] = impedanceMatrix[i][j] + 
                                      factorIJ * impedanceMatrix[i][imageJ] + 
                                      factorJI * impedanceMatrix[imageI][j] + 
                                      factorIJ * factorJI * impedanceMatrix[imageI][imageJ];
                }
            }
        }
        
        // Update the impedance matrix
        impedanceMatrix = tempMatrix;
        
        return true;
    } catch (...) {
        return false;
    }
}

// Factory function
__declspec(dllexport) std::shared_ptr<SymmetryBoundaryManager> createSymmetryBoundaryManager() {
    return std::make_shared<SymmetryBoundaryManager>();
}

// Utility functions
SymmetryType stringToSymmetryType(const std::string& typeStr) {
    std::map<std::string, SymmetryType> typeMap = {
        {"none", SymmetryType::NONE},
        {"symmetric", SymmetryType::SYMMETRIC},
        {"antisymmetric", SymmetryType::ANTISYMMETRIC},
        {"pec", SymmetryType::PEC},
        {"pmc", SymmetryType::PMC},
        {"asymmetric_excitation", SymmetryType::ASYMMETRIC_EXCITATION}
    };
    
    std::string lowerTypeStr = typeStr;
    std::transform(lowerTypeStr.begin(), lowerTypeStr.end(), lowerTypeStr.begin(), ::tolower);
    
    auto it = typeMap.find(lowerTypeStr);
    if (it != typeMap.end()) {
        return it->second;
    }
    
    return SymmetryType::NONE;  // Default
}

std::string symmetryTypeToString(SymmetryType type) {
    switch (type) {
        case SymmetryType::NONE:
            return "none";
        case SymmetryType::SYMMETRIC:
            return "symmetric";
        case SymmetryType::ANTISYMMETRIC:
            return "antisymmetric";
        case SymmetryType::PEC:
            return "pec";
        case SymmetryType::PMC:
            return "pmc";
        case SymmetryType::ASYMMETRIC_EXCITATION:
            return "asymmetric_excitation";
        default:
            return "none";
    }
}

PlaneOrientation stringToPlaneOrientation(const std::string& planeStr) {
    std::map<std::string, PlaneOrientation> planeMap = {
        {"xy", PlaneOrientation::XY_PLANE},
        {"yz", PlaneOrientation::YZ_PLANE},
        {"xz", PlaneOrientation::XZ_PLANE},
        {"z=0", PlaneOrientation::XY_PLANE},
        {"x=0", PlaneOrientation::YZ_PLANE},
        {"y=0", PlaneOrientation::XZ_PLANE}
    };
    
    std::string lowerPlaneStr = planeStr;
    std::transform(lowerPlaneStr.begin(), lowerPlaneStr.end(), lowerPlaneStr.begin(), ::tolower);
    
    auto it = planeMap.find(lowerPlaneStr);
    if (it != planeMap.end()) {
        return it->second;
    }
    
    return PlaneOrientation::XY_PLANE;  // Default
}

std::string planeOrientationToString(PlaneOrientation plane) {
    switch (plane) {
        case PlaneOrientation::XY_PLANE:
            return "XY_PLANE";
        case PlaneOrientation::YZ_PLANE:
            return "YZ_PLANE";
        case PlaneOrientation::XZ_PLANE:
            return "XZ_PLANE";
        default:
            return "XY_PLANE";
    }
}

} // namespace SymmetryBoundary